Efficient solutions of self-consistent mean field equations for dewetting and 

electrostatics in nonuniform liquids 



Zhonghan HrQ and John D. Weeks 
Institute for Physical Science and Technology and Department of Chemistry and Biochemistry, 
University of Maryland, College Park, Maryland 20742 
(Dated: September 17, 2010) 

We use a new configuration-based version of linear response theory to efficiently solve self- 
consistent mean field equations relating an effective single particle potential to the induced density. 
The versatility and accuracy of the method is illustrated by applications to dewetting of a hard 
sphere solute in a Lennard- Jones fluid, the interplay between local hydrogen bond structure and 
electrostatics for water confined between two hydrophobic walls, and to ion pairing in ionic solutions. 
Simulation time has been reduced by more than an order of magnitude over previous methods. 



Mean field theories have long provided a physically 
suggestive and qualitatively useful description of struc- 
ture, thermodynamics, and phase transitions in con- 
densed matter systems. In these approaches certain 
long-ranged components of the intermolccular interac- 
tions are replaced by an effective single particle potential 
or "molecular field" that depends self-consistently on the 
nonuniform density the field itself induces Recent 
work has shown that this method can produce exception- 
ally accurate results for models of simple and molecular 
liquids, ionic fluids, and water provided that i) all the in- 
termolecular interactions averaged over are slowly vary- 
ing at typical nearest neighbor distances where strong 
local forces dominate and ii) an accurate description of 
the density induced by the effective field is used 
When both these conditions are satisfied we call the re- 
sulting approach Local Molecular Field (LMF) theory. 

LMF theory can be viewed as a mapping that relates 
the structure and thermodynamics of the full long-ranged 
system to those of a simpler "mimic system" with trun- 
cated intermolecular interactions but in an effective or 
restructured field that accounts for the averaged effects of 
the long-ranged interactions @, Q . Analyzing the short- 
ranged mimic system may offer particular advantages 
for systems with Coulomb interactions, since standard 
particle-mesh Ewald sum treatments do not scale well 
in massively-parallel simulations [f|. But to realize the 
full potential of LMF theory as a practical tool in com- 
puter simulations and in qualitative analysis, we must ef- 
ficiently determine the effective field. We show here that 
when condition i) above is satisfied we can use a new 
configuration-based version of linear response theory to 
satisfy condition ii) as well. This leads to a highly ac- 
curate and efficient way to solve the self-consistent LMF 
equation. These ideas may also help solve self-consistent- 
field equations that appear in many other contexts. 

To illustrate the method in its simplest form, we first 
study dewetting of a single hard sphere "solute" particle 
in a Lennard- Jones (LJ) fluid. However LMF theory pro- 
vides a unified perspective and only minor changes are 
needed to describe more complex systems with Coulomb 
interactions. Here we focus on the coupling between 



dewetting and electrostatics in water confined by hy- 
drophobic walls and on ion pairing in ionic solutions. 
Long-ranged forces play a key but subtle role in all these 
examples and require a very accurate solution of the LMF 
equation. 

Consider a nonuniform LJ fluid in the presence of the 
field 4>q (r) generated by the hard sphere solute, as de- 
scribed below. The LJ pair potential is separated by 
the sign of the force into a short-ranged repulsive core 
and a longer-ranged perturbation part 0: ULj(rij) — 
uo(rij) + U\(rij). U\ contains only slowly- varying attrac- 
tive forces and thus is suitable for LMF averaging. The 
LJ mimic system has truncated intermolccular interac- 
tions uo(rij) in the presence of an effective or restruc- 
tured field </>r(t), chosen in principle so that the nonuni- 
form density pr,(t; [</>r]) = (p(r, R)) , in the mimic sys- 
tem (indicated by the subscript R) equals the density 
p(r; [4>o]) in the full system J2t5| ■ Here ( denotes a 
normalized ensemble average in the mimic system in the 
presence of </>r, R = {r^} denotes a microscopic config- 
uration of all particles, and p(r, R) = X^=i S(r — r.j) is 
the microscopic configurational density. 

As shown in detail elsewhere , and the asso- 
ciated equilibrium density PR,(r; [4>r]) can be accurately 
determined by solving the self-consistent LMF equation 

R (r) = Mr) + J dv' PK {v'- [fo]) Ul (|r - r'|) + C, (1) 

where C is a constant setting the zero of energy [§] . 

In earlier work 0, [If , the LMF equation was solved by 
straightforward iteration, using computer simulations to 
accurately determine the density induced by the given 
field at each iteration. However, even with a good initial 
estimate </>o for the effective field, determining the density 
induced by the remaining changes in the field 0ri = </>r — 
(j)Q from further iterations required more simulations. By 
rewriting the density in Eq. (TTJ) so that the dependence on 
0ri appears inside the ensemble average, we show here 
that the LMF equation can be iterated to self-consistency 
with no new simulations required. 

This is very easy to do formally. We refer to the system 
with effective field (f>o as a "trial" system and note that 
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the total potential energy associated with the correction 
</>ri in a configuration R is given by 
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Similarly defining J7q(R) as the total intermolecular po- 
tential energy in configuration R we have exactly 
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The idea behind Eq. usually with a configuration 
dependent function F(R) replacing p(r, R) and different 
choices of trial and full systems, is very well known and 
serves as the basis for perturbation and weighted his- 
togram methods for thermodynamic properties Q . Here 
we use it to determine changes of the density in the LMF 
equation as the iteration proceeds to self-consistency. 

Averaging over exponentials as in Eq. (|3|) is usu- 
ally problematic, and in most applications sophisticated 
techniques like umbrella sampling are needed to obtain 
good statistics 0]. However, condition i), required for 
the quantitative validity of LMF theory itself, ensures 
that only slowly varying intermolecular forces appear 
in <i> R1 (R,). This suggests it may not be difficult to 
find a trial field <f>o for which the exponential remain- 
der e -/3 * R1 ( R ' varies sufficiently slowly over most relevant 
configurations R that a simple direct average is accurate. 

A linearization of Eq. (j3|) provides a way to test for this 
condition. Defining <y$ R1 (R) = $ri(R) - (I»r,i(R))i 

and Sp(r, R) = p(r, R) - (p(r, R)) i , we find 



<p(r,R)> ~ <p(r,R)), -0(*p(r,R)**Ri(R) 



(4) 

This is equivalent to the usual linear response formula 

-P j dr' (Sp(r,R)Sp(v',K)) $ J m (v') (5) 

on using Eq. @ . However the linear response (LR) func- 
tion \6p(r, R)5p(r', R))^ depends on r and r' separately 
in a nonuniform system and is usually too complicated to 
be determined directly. In contrast, the averages in Eqs. 
dD) and ([3]) simply reweight the bin histograms used to 
determine the nonuniform density and depend on r alone. 
They can be easily calculated using the saved configura- 
tions of a single well-chosen trial system. 

For very slowly varying <I>ri(R), the system is in a lin- 
ear regime where essentially identical structural changes 




FIG. 1. LMF/LR/EXP treatment of the drying transition 
induced by a hard sphere solute in a LJ fluid for T = 0.85, 
p B — 0.70. (a) The RDF around a hard sphere cavity with 
R c = 2 for the full LJ fluid, the SCA trial system, and the 
converged LMF system, (b) Trial fields in the linear regime 
(solid and dot-dashed curves) and the final self-consistent field 
(circles and solid line) 



arise from the exponential (EXP) form in Eq. ([3]) or the 
LR form in Eq. By requiring that a trial system 

gives the same result on iterating the LMF equation us- 
ing both forms, we have a conservative and objective test 
for an accurate solution. We refer to this as the LR-LMF 
method, since in general Eq. (|4]) proves most useful. To 
our knowledge, the advantages of the configuration-based 
LR formula in Eq. (j4j have not been exploited before. 

We first study a hard sphere solute in a LJ fluid at a 
state near the triple point [l(| . The solute is represented 
by an external field 4>o{r) that is infinite inside a cavity 
of radius R c and zero otherwise. In the simplest "strong 
coupling approximation" (SCA) to LMF theory, often 
successfully used in perturbation theories of uniform liq- 
uids 0, all effects of the slowly varying u\ on the molec- 
ular structure are ignored and the effective field 4>r(t) is 
approximated by the bare hard sphere field (f>o(r). When 
R c (measured in units of olj) is unity, attractive forces 
nearly cancel and the SCA gives good results, correctly 
predicting an oscillatory density response with a large 
density maximum at contact. 

However, as the cavity radius increases, particles near 
the cavity experience unbalanced attractive forces from 
LJ particles further away that reduce the contact den- 
sity. Strong reduction already is seen for R c = 2. 
The circles in Fig. la give results of computer simu- 
lations for the radial distribution function (RDF) 
g(r; [cj>o]) = p(r; [4>o\)l P B induced by the hard sphere so- 
lute in the full LJ fluid with bulk density p B . This ex- 
hibits an essentially structureless profile with a contact 
value of unity. The solid curve in Fig. la gives <?R,(r; [</>r]), 
the RDF in the mimic system. The excellent agreement 
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between the LMF prediction and the full LJ g(r; [4>o]) 
shows that LMF theory quantitatively captures the dry- 
ing effect. This contrasts with the dashed line in Fig. 
la, the density in the truncated system induced by the 
bare field 4>q. The failure of this SCA prediction illus- 
trates the general need in most nonuniform systems for 
a proper self-consistent solution of the LMF equation. 

The self-consistent field </>r(t") was obtained from solv- 
ing Eq. (fTJ) using two different trial fields shown in Fig. 
lb. Both trial systems are in the linear regime and give 
the same final result. The solid curve results from a pre- 
vious trial simulation based on the SCA that used the 
bare 4>o( r )- As Fig. la suggests, this represents a very 
poor initial guess, and the EXP form ([3]) produced very 
noisy data, indicating poor overlap between the SCA and 
final LMF configurations. However Eq. ((4]) gave much 
smoother data and its use in the LMF equation gave the 
solid curve in Fig. lb as the output field. Using this as 
a second trial field generates a converged solution of the 
LMF equation using either Eqs. (01 or 

Since U\ is slowly-varying, using relatively crude ap- 
proximations to the asymmetric density in the LMF 
equation can often give a better trial field than the sim- 
ple SCA. Thus in Fig. lb, the dot-dashed trial field was 
found by approximating the density by a step function 
that vanishes inside the cavity and equals/? 5 outside. 

LMF theory proves even more useful [J, [5j when ap- 
plied to systems with Coulomb interactions in the pres- 
ence of an electrostatic potential V(r) arising from a fixed 
external charge distribution p q xt (r'). The basic Coulomb 
interaction 1/r = Vq(t) + vi(r) is separated into short- 
and long-ranged components, where v\{r) is proportional 
to the electrostatic potential arising from a normalized 
Gaussian charge distribution with width a, 
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An advantage of the Coulomb separation is that a can be 
chosen specifically in different applications so that con- 
dition i) is very well satisfied. 

When all charges in the system (both fixed and mobile) 
are separated using the same er and all other intermolec- 
ular interactions remain unchanged, LMF theory then 
gives a mapping to a Coulomb mimic system where all 
1/r interactions are replaced by the short-ranged vo(r) 
and there is a restructured electrostatic potential Vr that 
satisfies the Coulomb LMF equation 



Vft(r) = V (r) + 



dr 'PR.tot 



(r';[V R ]W(|r-r'|)+C. (7) 



Here Vb(r) is the short-ranged part of the external po- 
tential, given by the convolution of Vo(r) with the fixed 
charge density, and tot (r'; [Vr]) is the total equilib- 
rium charge density from both fixed and mobile charges. 
An alternate form of Eq. ([7]) better relates LMF the- 




FIG. 2. Charge densities for water confined by hydrophobic 
LJ walls pl| centered at z = ±2.25 nm for the SCA (dashed 
line), mimic system (solid line) and full system (circles). The 
inset shows the corresponding polarization potentials [ill ]. 



ory to conventional electrostatics Noting the convo- 
lution defining v\, we see that the slowly- varying part 
VRi(r) = Vr(i") — Vo(r) of the restructured potential 
in ([7]) exactly satisfies Poisson's equation but with a 
Gaussian-smoothed charge density j0p t ° tot , given by the 
convolution of p Rtot with the Gaussian in Eq. ©. 

We now apply the Coulomb LMF Eq. ([7|) to the ex- 
tended simple point charge (SPC/E) model for water fl2| . 
As suggested by the SCA, we first consider a Gaussian- 
truncated (GT) model, where all the 1/r interactions 
from charges in SPC/E water are replaced by v (r) and 
we ignore all structural effects from V\. Bulk GT wa- 
ter with a = 0.45 nm gives excellent results for both 
atom-atom and dipolc-angle correlation functions when 
compared to the full SPC/E model, with Coulomb inter- 
actions treated by Ewald sums 

Moreover, as Fig. 2 shows, when GT water is confined 
between two hydrophobic LJ walls as defined in 1131. the 
charge density Pq(z) determined by simulations [ll| us- 
ing bare wall fields with Vo = seems to capture most 
qualitative features of the dipole layer, a characteristic 
property of water near extended hydrophobic interfaces 
[5j, [l3j . Only small differences in the peak heights are 
visible when compared with full SPC/E water, simulated 
using the slab-corrected Ewald 3D method 11-411 . 

Nevertheless, as has long been recognized [l5lll6j|. ma- 
jor errors are seen in the polarization potential felt by a 
test charge, given by integrating Poisson's equation: 



$ P oi(z) 



dz' 



L/2 



L/2 



dz"p q {z"). 



(8) 



As shown in the inset of Fig. 2 for the full SPC/E model, 
$ P oi should reach a plateau in the central bulk region, 
and GT water fails dramatically in this respect. 

A self-consistent solution of Eq. ([7]) yields a very accu- 
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excellent results have also been obtained for ion pairing 
in ionic solution models using the LR-LMF method 11 



-1 1 
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FIG. 3. Gaussian-smoothed charge densities with a = 0.45 
nm. Labels are the same as in Fig. [2] The inset gives the 
effective field Va(r). Note that the bare V(r) = in this case. 



rate charge density that corrects all such failures. The in- 
set in Fig. 3 gives the converged Vr(z). Gaussian smooth- 
ing of the bare charge density as dictated by LMF the- 
ory averages over the simulation noise and local structure 
and reveals the much smaller coherent long- wavelength 
features that control the electrostatics [B| . The smoothed 
charge density quickly decays to a neutral bulk in Fig. 3 
for both SPC/E water and LMF theory, while a true bulk 
never forms in GT water, causing the very poor polar- 
ization potential in Fig. 2. 

These results highlight the advantages of the new LR- 
LMF method. It reduced the simulation time by more 
than an order of magnitude compared to use of the stan- 
dard iteration method [j| or to use of the slab-corrected 



Ewald 3D method [14j. When using the standard iter- 



ation method, it proved very difficult to distinguish be- 
tween equilibrium charge density fluctuations, present for 
any given field, and the desired changes in the charge 
density as the iteration proceeded to self-consistency. To 
obtain convergence, earlier workers had to use a large 
value of a = 0.6 nm and the charge density at each iter- 
ation was taken as the average over a set of 10 parallel 
simulations with different initial conditions 

In contrast, the results in Figs. 2 and 3 were obtained 
from only two very short trial simulations [111 ], starting 
first with the SCA Vo = 0, and using the smaller bulk 
a = 4.5 nm. As before, the EXP form ([3]) was very noisy 
when used with SCA configurations, but the LR form 
((4]) was much better behaved. It generated a second trial 
field in the linear regime that was virtually identical to 
the final self-consistent field shown in the inset of Fig. (3). 
The remaining effects of equilibrium fluctuations in the 
mimic system show up as small variations (about the size 
of the circle symbol) in the bulk value of the polarization 
potential in Fig. 2 or the heights of the charge density 
peaks in Fig. 3. As shown in the supplementary material, 



Reaction field (RF) truncations of Coulomb interac- 
tions have recently been used in massively-parallel simu- 
lations of biological systems to permit much faster sim- 
ulations Q. LR-LMF theory may provide a promising 
linear-scaling alternative in which the effective field cor- 
rects known problems like those illustrated in Fig. 2 aris- 
ing from simple RF or SCA truncations. More generally, 
we believe the efficient solutions generated by the LR- 
LMF method establish the full power of the basic mean 
field picture for both quantitative and qualitative anal- 
ysis of a wide range of electrostatic and dielectric phe- 
nomena in nonuniform liquids. 
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Jocelyn Rodgers for very helpful remarks. 
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L SIMULATION DETAILS FOR DEWETTING OF THE LJ FLUID AROUND A HARD SPHERE SOLUTE AND FOR WATER 

CONFINED BETWEEN HYDROPHOBIC WALLS 

In the following, equation and figure numbers refer to the main paper, unless otherwise indicated. Our Monte Carlo (MC) 
simulation for N = 2363 LJ particles with a hard sphere solute at a state near the triple point with T = 0.85 and p B = 0.70 
follows the work of Huang and Chandler [Q]]. Configurations of a short trial simulation of 10000 MC steps based on the SCA 
and Eq. (4) were used in the LMF Eq. (1) to generate the trial field given by the solid line in Fig. (lb). A second trial simulation 
of 20000 steps using that field generated the final self-consistent LMF. Both runs were much shorter than the final production 
run of 200000 steps needed to determine the smooth profile of the full LJ system given by the circles in Fig. (la). 

Our molecular dynamics (MD) simulation of SPC/E water confined between hydrophobic LJ walls used the same system 
parameters and a modified DLPOLY2. 1 8 MD package Q2Q similar to that described in detail in [3|], with further modifications to 
save configurations needed for the LR-LMF theory. Two trial simulations with a = 4.5 nm each of 1 ns length and starting from 
the SCA were used to generate the data in Figs. (2) and (3). This was more than an order of magnitude faster than analogous 
results obtained using the standard iteration method |@]. Moreover, as discussed in detail in J3l, a large value of a = 0.6 nm had 
to be used in the standard iteration method, because of problems in distinguishing equilibrium charge density fluctuations with 
a given estimate for the effective field from the desired iterative changes in the charge density. Even after averaging over results 
of 10 parallel simulations for each iteration, convergence was not found with a = 4.5 nm. 

In contrast, the LR-LMF method uses the same set of equilibrium trial system configurations to determine the remaining 
changes in the charge density as self-consistency is achieved by iterating the LMF equation. Convergence was found using the 
same a = 4.5 nm used for the bulk, and required only two short trial simulations. 

H. CASE STUDY OF DILUTE IONIC SYSTEM 

To further illustrate the versatility of the LR-LMF method we discuss here results for a dilute bulk ionic system at moderately 
strong coupling. This presents a challenge to theory, since the coupling is strong enough that traditional methods based on 
Poisson-Boltzmann theory will fail, but weak enough that the SCA alone will not be accurate unless a very large a is used. 
Indeed, previous applications of LMF theory to bulk ionic solution models considered only strong-coupling states and used such 
a large a that the SCA alone was accurate yfl. The present calculation represents the first self-consistent solution of the full 
LMF equation for a general ionic solution model in this difficult regime. As we will see only a single simulation in this case is 
needed to achieve full self-consistency. 

The system is a size-asymmetric primitive model consisting of charged LJ molecules. The LJ interaction parameters are 

(j + = 4.0A, cr_ = 3.0A, o"-| = 3.5A, e + = O.lkcal/mol, e_ = 0.4kcal/mol and = 0.2kcal/mol. The cation and anion 

has unit charge and mass of 12 and 30 respectively. The cutoff distance for both the LJ and Ewald short ranged interaction in the 
system is 15 A. The system has 1024 pairs of ions in a cubic box of size 62 A with T = 5000 K. The MD simulation of the full 
system using Ewald sums used 90000 steps. 
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FIG. 1 . The restructured field with a cation fixed at the origin 

Results of MD simulations for the full system and for the LR-LMF method using the SCA as the trial system are presented 
here. Fig. (1) shows the restructured field with a cation fixed at the origin. The charge density around the cation is shown 
in Fig. (2). LR-LMF and EXP-LMF results overlap starting from SCA configurations and we show only the former in the 
graphs. Significant ion pairing is indicated by the strong first peak in the cation-anion correlation function in Fig. (3). The SCA 
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FIG. 2. The charge density around a cation 
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FIG. 3. The cation-anion correlation function 



simulation used a = 3.5A, a value for which the SCA correlation function had noticeable errors, as seen in Fig. (3). However 
the trial system based on SCA configurations is in the linear regime, and no further simulations are needed to get the converged 
results shown here from iteration of the LMF equation. 

As explained in |6L it is useful to rewrite the LMF equation in a more numerically stable form for iteration so that small errors 
in the decay of the charge density are not amplified when convoluted with the long-ranged v\ . We used the specific method 
described in Appendix E of 0]. 

It is gratifying that LR-LMF theory allows for such accurate results using a small value of a, since this permits fast mimic 
system simulations. However, care must be taken, since the underlying LMF theory itself will fail if too small a a is used, and 
the small differences still visible for the LR-LMF and full correlation functions suggests we are close to the minimum accurate 
value. 
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